1 Loading Packages + Data

load(here("jk_code", "JK_cleanMD.rds"))

2 Removing mt and Rp genes

SO4 <- subset(SO4, features = grep("^mt-|^Rp", rownames(SO4), invert = TRUE, value = TRUE))

3 Reclustering and finding Markers

SO4 <- SCTransform(SO4) %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:30) %>%
    FindClusters(resolution = 0.3) %>%
    RunUMAP(dims = 1:30)
## Running SCTransform on assay: RNA
## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 22376 by 11431
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## There are 2 estimated thetas smaller than 1e-07 - will be set to 1e-07
## Found 306 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 22376 genes
## Computing corrected count matrix for 22376 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 20.65131 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Set default assay to SCT
## PC_ 1 
## Positive:  Umod, Egf, Tmem52b, Sult1d1, Sostdc1, Klk1, Prdx5, Fabp3, Cox6c, Wfdc2 
##     Foxq1, Ly6a, Wfdc15b, Slc25a5, Cox5b, Atp5g1, Krt7, Ckb, Cox4i1, Atp1a1 
##     Cldn19, Ndufa4, Cox8a, Ggt1, Atp5b, Chchd10, Cox6b1, Atp1b1, Uqcr11, Cox7a2 
## Negative:  Pappa2, Malat1, Gm37376, Zfand5, Mir6236, Robo2, Aard, Neat1, Wwc2, Pde10a 
##     Itga4, Nadk2, Col4a4, Nos1, Col4a3, Sgms2, S100g, Irx1, Ptgs2, Etnk1 
##     Ramp3, Sdc4, Srrm2, Itprid2, Hnrnpa2b1, Tmem158, Bmp3, Ranbp3l, Camk2d, Zbtb20 
## PC_ 2 
## Positive:  Egf, Mir6236, Malat1, Gm37376, Umod, CT010467.1, Tmem52b, Nme7, Neat1, Slc12a1 
##     Gm24447, Kcnq1ot1, Atp1b1, Etnk1, Wnk1, Gm28438, Lars2, Sult1d1, Dst, Foxq1 
##     Slc5a3, Atrx, Pnisr, Syne2, Ivns1abp, Hnrnpa2b1, Gls, Srrm2, Gm23935, Col4a4 
## Negative:  Gm22133, Ubb, Fth1, Mt1, Ftl1, Hspa1a, Ldhb, H3f3b, Eif1, Car15 
##     Prdx1, Cd63, Hspa1b, Fos, Junb, Cd9, Mt2, Fxyd2, Mpc2, Mgst1 
##     Clu, Tmem213, Bsg, Jun, Jund, Actb, Spp1, Aldoa, Itm2b, Btg2 
## PC_ 3 
## Positive:  Fos, Junb, Jun, Hspa1b, Btg2, Hspa1a, Egr1, Zfp36, Atf3, Ier2 
##     Fosb, Klf6, Socs3, Klf2, Jund, Dnajb1, Dusp1, Gadd45g, Rhob, Tsc22d1 
##     Gadd45b, Ubc, Malat1, Cebpd, H1f2, Klf4, Gm42793, H2bc4, H3f3b, Nr4a1 
## Negative:  Pappa2, Car15, Gm22133, Cd9, Mgst1, Ldhb, Fth1, Aard, Itm2b, Mpc2 
##     Cd63, Tmem59, Bsg, Prdx1, Clu, Tmem213, Ramp3, Tmem158, Ftl1, Tmbim6 
##     Mcub, Sfrp1, Epcam, Ndfip1, Fxyd2, Tmem176a, Irx1, Ly6a, Mdh1, Rabac1 
## PC_ 4 
## Positive:  Pappa2, Aard, Umod, Tmem52b, Egf, Sult1d1, Tmem158, Mcub, Car15, Ptgs2 
##     Wwc2, Foxq1, Hsp90b1, Dctd, Iyd, Cd9, Ptprz1, Wnk1, Ramp3, Fth1 
##     Cdkn1c, Clu, Pth1r, Gm44120, Hspa5, Sec14l1, Ctsc, Calr, Col4a3, 1700028P14Rik 
## Negative:  Mt1, Gm8420, Apoe, Mt2, Gm8885, Gm28438, Aebp1, Gm24447, Sostdc1, Gpx6 
##     Gm28437, Neat1, Gm11808, Gm28037, Gm3511, Fgf9, Atp5md, Egfl6, Fxyd2, Gm45895 
##     Ivns1abp, Defb1, Gm11353, Atp5k, Gm12338, Gm10221, Cox6c, Car4, Gm6136, Ptger3 
## PC_ 5 
## Positive:  S100g, Gm8420, Gm8885, Actb, Tmem52b, Gm11808, Gm11353, Ppia, Tmsb10, Gm6136 
##     Gm3511, Gm12338, Gm9616, Gm10029, Igfbp7, Gm7206, Gm9843, Gm28037, Gm10221, Serf2 
##     Atp5md, Gm12254, Atp5e, Aard, Ndufb1-ps, Cox6c, Gm4332, Cebpd, Pappa2, Gm5526 
## Negative:  Hspa1a, Klk1, Hspa1b, Gm22133, CT010467.1, Car15, Gm13340, Fth1, Ptger3, Itm2b 
##     Atp1a1, Cldn10, Hspa8, Apoe, Cd63, Aebp1, Ftl1, Mt2, Gpx6, Uba52-ps 
##     Wfdc2, Jun, Fau, Mt1, Pik3r1, Lpl, Gm6061, Clcnkb, Sostdc1, Gm3362
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11431
## Number of edges: 377988
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8604
## Number of communities: 6
## Elapsed time: 1 seconds
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 14:48:16 UMAP embedding parameters a = 0.9922 b = 1.112
## 14:48:16 Read 11431 rows and found 30 numeric columns
## 14:48:16 Using Annoy for neighbor search, n_neighbors = 30
## 14:48:16 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 14:48:16 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpvLVyR4/file16105c2114f0
## 14:48:16 Searching Annoy index using 1 thread, search_k = 3000
## 14:48:17 Annoy recall = 100%
## 14:48:18 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 14:48:18 Initializing from normalized Laplacian + noise (using RSpectra)
## 14:48:18 Commencing optimization for 200 epochs, with 506010 positive edges
## 14:48:18 Using rng type: pcg
## 14:48:21 Optimization finished
SO.markers <- FindAllMarkers(SO4, only.pos = TRUE)
## Calculating cluster 0
## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
SO.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)
SO.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10
DoHeatmap(SO4, features = top10$gene) + NoLegend()
## Warning in DoHeatmap(SO4, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay: Phkg1,
## Gm26982, Rtp4, Rsad2, Oasl2, Oasl1, Il3ra, Gm3500, Gm8257

DimPlot(SO4,)

DimPlot(SO4, split.by = "sample")

markers.to.plot2 <- c("Lrp2",         # PT
                      "Slc5a12",      # PT-S1
                      "Slc13a3",      # PT-S2
                      "Slc16a9",      # PT-S3
                      "Havcr1",       # Injured PT
                      "Epha7",        # dTL
                      "Cryab",        # dTL
                      "Cdh13",        # dTL1
                      "Slc14a2",      # dTL2
                      "Slc12a1",      # TAL
                      "Umod",         # TAL, DCT1
                      "Egf",          # TAL, DCT1,
                      "Cldn10",       # TAL
                      "Cldn16",       # TAL
                      "Nos1",         # MD
                      "Slc12a3",      # DCT
                      "Pvalb",        # DCT1
                      "Slc8a1",       # DCT2, CNT
                      "Aqp2",         # PC
                      "Slc4a1",       # IC-A
                      "Slc26a4",      # IC-B
                      "Nphs1",        # Podo
                      "Ncam1",        # PEC
                      "Flt1",         # Endo
                      "Emcn",         # Glom Endo
                      "Kdr",          # Capillary Endo
                      "Pdgfrb",       # Perivascular
                      "Pdgfra",       # Fib
                      "Piezo2",       # Mesangial
                      "Acta2",        # Mural
                      "Ptprc",        # Immune
                      "Cd74",         # Macrophage
                      "Skap1",        # B/T Cells 
                      "Upk1b",        # Uro
                      "Top2a",         # Proliferation
                      "Cldn5",
                      "Jun",
                      "Fosb"
)
                      
DotPlot(SO4,
features = markers.to.plot2,
dot.scale = 8,
dot.min = 0,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()
## Warning: The following requested variables were not found: Slc5a12, Slc13a3,
## Havcr1, Kdr, Pdgfrb, Pdgfra, Piezo2, Upk1b

markers.to.plot1 <- c("Ramp3",         # PT
                      "Egf",      # PT-S1
                      "Fos",
                      "Cxcl10")      # PT-S3
                
DotPlot(SO4,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()

DimPlot(SO4)

DimPlot(SO4,group.by = "treatment",split.by = "treatment")

DimPlot(SO4,split.by = "sample")

4 Viewing the

markerstoplot <- c("Hspa8", "Ptger3", "Gxp4","Ckb","Atp4a","S100g")

DotPlot(SO4,
features = markerstoplot,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()
## Warning: The following requested variables were not found: Gxp4

DimPlot(SO4,split.by = "sample")

DimPlot(SO4,split.by = "sample",group.by = "seurat_clusters")

5 Subclass

  SO4@meta.data <- SO4@meta.data %>% 
  mutate(subclass_MD = dplyr::case_when(
    seurat_clusters == 0  ~ "type_1",
    seurat_clusters == 1  ~ "type_1",
    seurat_clusters == 2  ~ "type_2",
    seurat_clusters == 3  ~ "type_1",
    seurat_clusters == 4  ~ "type_3",
    seurat_clusters == 5  ~ "type_4",

  ))

SO4@meta.data$subclass_MD <- factor(SO4@meta.data$subclass_MD , levels = c("type_1", "type_2", "type_3", "type_4"))

Idents(SO4) <- SO4@meta.data$subclass_MD

DimPlot(object = SO4, reduction = "umap", group.by = "subclass_MD", label = TRUE)

DimPlot(object = SO4, reduction = "umap", label = TRUE)

Idents(SO4) <- "subclass_MD"

DimPlot(SO4,split.by = "sample",group.by = "seurat_clusters")

6 Viewing the DEGs between each type.

DotPlot(SO4,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()

7 Saving File

save(SO4, file = here("jk_code", "JK_remove_mtrp.rds"))